program test_mesh

  use flogger
  use fiona
  use string
  use const_mod
  use cubed_sphere_mesh_mod
  use cubed_sphere_array_mod

  implicit none

  integer, parameter :: nc = 48
  type(cubed_sphere_mesh_type) mesh
  type(cubed_sphere_array_type) q
  integer is, ie, js, je, start(3), count(3)
  character(5) :: cell_dims (3) = ['x    ','y    ','panel']
  character(5) :: vtx_dims  (3) = ['xi   ','yi   ','panel']
  character(5) :: xedge_dims(3) = ['xi   ','y    ','panel']
  character(5) :: yedge_dims(3) = ['x    ','yi   ','panel']
  real(r8) :: sum_area = 0

  integer ipn

  call mesh%init(proj_type=1, nc=nc, nz=1, hw=2, ipn=1, ids=1, jds=1, ndx=5, ndy=5)
  call q%init(mesh, 'q', '', '', loc=grid_loc_cell)

  do ipn = 1, 6
    if (mesh%panels(ipn)%active) then
      associate (q => q%panels(ipn)%values)
      q(:,:,:) = ipn
      end associate
    end if
  end do

  call fiona_init()

  call fiona_create_dataset('f', desc='Cubed-sphere mesh', file_path='c' // to_str(nc) // '.nc')
  call fiona_add_dim('f', 'panel', size=6   )
  call fiona_add_dim('f', 'x'    , size=nc  )
  call fiona_add_dim('f', 'y'    , size=nc  )
  call fiona_add_dim('f', 'xi'   , size=nc+1)
  call fiona_add_dim('f', 'yi'   , size=nc+1)
  call fiona_add_dim('f', 'four' , size=4   )
  call fiona_add_var('f', 'cell_lon'    , long_name='Longitude of cells'      , units='degree_east' , dim_names=cell_dims, dtype='r8', missing_value=r8_missing_value)
  call fiona_add_var('f', 'cell_lat'    , long_name='Latitude of cells'       , units='degree_north', dim_names=cell_dims, dtype='r8', missing_value=r8_missing_value)
  call fiona_add_var('f', 'vtx_lon'     , long_name='Longitude of vertices'   , units='degree_east' , dim_names=vtx_dims , dtype='r8', missing_value=r8_missing_value)
  call fiona_add_var('f', 'vtx_lat'     , long_name='Latitude of vertices'    , units='degree_north', dim_names=vtx_dims , dtype='r8', missing_value=r8_missing_value)
  call fiona_add_var('f', 'cell_area'   , long_name='Area of cells'           , units='m2'          , dim_names=cell_dims, dtype='r8', missing_value=r8_missing_value)
  call fiona_add_var('f', 'cell_subarea', long_name='Sub-area of cells'       , units='m2'          , dim_names=['four ','x    ','y    ', 'panel'], dtype='r8', missing_value=r8_missing_value)
  call fiona_add_var('f', 'cell_lx'     , long_name='Cell length along x axis', units='m'           , dim_names=cell_dims, dtype='r8', missing_value=r8_missing_value)
  call fiona_add_var('f', 'cell_ly'     , long_name='Cell length along y axis', units='m'           , dim_names=cell_dims, dtype='r8', missing_value=r8_missing_value)
  call fiona_add_var('f', q%name, long_name=q%long_name, units=q%units, dim_names=cell_dims, dtype='r8')
  call fiona_start_output('f', 0.0d0, new_file=.true.)
  do ipn = 1, 6
    if (mesh%panels(ipn)%initialized) then
      is = mesh%panels(ipn)%full_ids; ie = mesh%panels(ipn)%full_ide
      js = mesh%panels(ipn)%full_jds; je = mesh%panels(ipn)%full_jde
      start = [is, js, ipn]
      count = [mesh%panels(ipn)%full_ndx, mesh%panels(ipn)%full_ndy, 1]
      call fiona_output('f', 'cell_lon'    , mesh%panels(ipn)%cell%lon_deg(  is:ie,js:je), start=start, count=count)
      call fiona_output('f', 'cell_lat'    , mesh%panels(ipn)%cell%lat_deg(  is:ie,js:je), start=start, count=count)
      call fiona_output('f', 'cell_area'   , mesh%panels(ipn)%cell%area   (  is:ie,js:je), start=start, count=count)
      call fiona_output('f', 'cell_lx'     , mesh%panels(ipn)%cell%lx     (  is:ie,js:je), start=start, count=count)
      call fiona_output('f', 'cell_ly'     , mesh%panels(ipn)%cell%ly     (  is:ie,js:je), start=start, count=count)
      call fiona_output('f', 'cell_subarea', mesh%panels(ipn)%cell%subarea(:,is:ie,js:je), start=[1,is,js,ipn], count=[4,mesh%panels(ipn)%full_ndx,mesh%panels(ipn)%full_ndy,1])
      sum_area = sum_area + sum(mesh%panels(ipn)%cell%area(is:ie,js:je))
      is = mesh%panels(ipn)%half_ids; ie = mesh%panels(ipn)%half_ide
      js = mesh%panels(ipn)%half_jds; je = mesh%panels(ipn)%half_jde
      start = [is, js, ipn]
      count = [mesh%panels(ipn)%half_ndx, mesh%panels(ipn)%half_ndy, 1]
      call fiona_output('f', 'vtx_lon' , mesh%panels(ipn)%vtx %lon_deg(is:ie,js:je), start=start, count=count)
      call fiona_output('f', 'vtx_lat' , mesh%panels(ipn)%vtx %lat_deg(is:ie,js:je), start=start, count=count)

      is = mesh%panels(ipn)%full_ids; ie = mesh%panels(ipn)%full_ide
      js = mesh%panels(ipn)%full_jds; je = mesh%panels(ipn)%full_jde
      start = [is, js, ipn]
      count = [mesh%panels(ipn)%full_ndx, mesh%panels(ipn)%full_ndy, 1]
      call fiona_output('f', q%name, q%panels(ipn)%values(is:ie,js:je,1), start=start, count=count)
    end if
  end do
  call fiona_end_output('f')

  ! if (abs((sum_area - 4 * pi * radius**2) / sum_area) > 1.0e-12) then
  !   call log_error('Cell area is not matched with spherical area!')
  ! end if

  call fiona_final()

end program test_mesh
